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Inversion of marine heat flow measurements by expansion of the 
temperature decay function 
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SUMMARY 

Marine heat flow data, obtained with a Lister-type probe, consists of two temperature 
decay curves, frictional and heat pulse decay. Both follow the same physical model of a 
cooling cylinder. The mathematical model describing the decays is nonlinear as to the 
thermal sediment parameters thus a direct inversion is not possible. To overcome this 
difficulty, the model equations are expanded using a first order Taylor series. The lin- 
earised model equations are used in an iterative scheme to invert the temperature decay 
for undisturbed temperature and thermal conductivity of the sediment. The inversion 
scheme is tested first for its theoretical limitations using synthetic data. Inversion of 
heat flow measurements obtained during a cruise of R/V SONNE in 1996 and needle 
probe measurements in material of known thermal conductivity show that the algorithm 
is robust and gives reliable results. The program can be obtained from the authors. 
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Introduction 

The growing interest in global energy and geochemical 
fluxes in the oceans and the availability of multi-penetration 
heat flow probes have has led to increased attention in ma- 
rine geothermal studies, which are mainly focused on the 
sediment covered areas of ridges and ridge flanks. Studies 
of accretionary wedges that incorporate heat flow investi- 
gations also reflect the need for additional constraints to 
model and understand the dewatering process of accreted 
sediments. This increased interest in marine heat flow data 
has helped to improve the existing measurement technique 
in two ways: The violin bow type h eat probe instrument, as 
described in Hyndman e t alj d!979h . has evolved over two 
decades of intensive use to a mature, mechanically robust 
instrument which now can be used in a routine way. Rapid 
electronic development led to an increased temperature res- 
olution of 1 mK and allowed a larger number of sensors 
to be mounted on one suing due to increased storage ca- 
pacity. Both developments now permit multi-penetration de- 
ployments (measurements in a pogo-style fashion) of up to 
24 hours per station. Advances in marine navigation now 
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Figure 1. Cross section of the sensor tube, housing the temperature sensors 
(thermistors) and the heating wires for in situ thermal conductivity mea- 
surements (see also Table[T}. 



permit very detailed studies of regional processes. World- 
wide coverage of the Global Positioning System (GPS) and 
Differential GPS help enormously in station keeping during 
a measurement. High accuracy in positioning of the probe 
is achieved b y using Long- or Short-Baseline navigation 
(UonesUl999l) . 

Figure Q] shows a simplified cross section of a sensor 
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Material 


Sensor string 


1 


0.2 


2.14 


copper, oil, glass 


Heating wire 


2 


13.2 


3.50 


Ni-Chrome 


Oil filling 


3 


0.2 


1.60 


oil, plastic 


steel tube 


4 


37.5 


3.9 


steel 



Ta ble 1. Thermal param eters of the sensor tube of figure [T] (af- 
ter lNagihara & Listed fl993l) ). 



string used by the Heat Flow Group of the University of 
Bremen (Germany) and the Pacific Geoscience Centre (Sid- 
ney, B.C., Canada). The sensors for temperature measure- 
ments in the sediment are housed in a hydraulic steel tube 
with an outer diameter of 8 mm and an internal diameter 
of 5 mm. Therma l param eters are shown in Table Q] (after 
Nagihara & Listed 0993)). The sensor string itself contains 
heating wires for in situ thermal conductivity measurements 
and wires leading to the thermistors for temperature mea- 
surements. The tube is oil-filled to improve the thermal con- 
tact between the temperature sensors and the tube. 

In order to illustrate the steps involved to process a heat 
flow measurement, Fig.[2]shows a typical data set. Measure- 
ments of undisturbed sediment tempera ture and cond uctivity 
follow the pulsed needle probe method flListeil IT 970h . When 
the sensor string penetrates the sediment the friction between 
sensor tube and sediment creates heat resulting in a temper- 
ature rise. The following temperature decay is recorded at 
10 s sample rate for a preset time span (7 minutes), after 
which a calibrated heat pulse of 20 s length is fired. The 
heat pulse decay is monitored for at least 7 minutes until the 
probe is pulled out of the sediment and the ship moves to the 
next measurement position with the probe about 200-300 m 
above the seafloor (pogo-style measurements). 

The processing of the raw measurements requires three 
steps: (1) determine undisturbed sediment temperatures 
from frictional decay (2) correct heat pulse decay for the re- 
maining effect of the frictional decay and (3) calculate ther- 
mal conductivities from heat pulse decay. The basic design 
of the processing of h eat flow measurements is outlined in 
iHvndman et al.l d 19791) wh ich was basic ally a manual proc e- 
dure based on the work of lListerl (1 19701) and lListerl ([1979). 

The increasing number of measurements per cruise in 
the past decade required a processing scheme which could 
easily be implemented on a personal computer for au- 
tomated processing of a large number of measurements. 
IVillinger & Davisl d 1987b published a pragmatic scheme 
(called HFRED), that minimises the misfit between mea- 
sured and model data in a least-squares sense by varying 
the effective origin time of penetration. Tests on numerically 
modelled data (synthetic measurement) with known param- 
eters showed that HFRED produced reliable and accurate re- 
sults. However, the scheme has two major deficiencies: 

(i) The thermal diffusivity used for the sediment is com- 
puted from t hermal conductivity ac cording to a relationship 
proposed by Hvndman et al.l d 1 979b . This relationship has 
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Figure 2. Example of a heat flow measurement. Only three out of 1 1 sensors 
are shown. 



never been validated by experimental data and will certainly 
vary with sediment type. 

(ii) The algorithm implemented in HFRED does not al- 
low analysis of errors of the calculated undisturbed sediment 
temperatures and in-situ thermal conductivities in a rigorous 
way; errors calculated by HFRED are always unrealistically 
low, co mpared to error estimates of about 5% from other 
studies (IHvndman et al.L 1 1 979t iListeii 1 1 970b . 

To overcome these deficiencies and to incorporate platform 
independent plotting routines, a mathematically sound in- 
version scheme of observed temperature decays was imple- 
mented using Matlab®, a widely used software package 
for numerical analysis. This allows creation of very com- 
pact code for the inversion, on-screen graphics and platform- 
independent plots. In addition, automated processing or re- 
processing of a large number of individual measurements 
is possible. The inversion of the integral describing the de- 
cay of a temperature pulse (see equation Q]) allows use of 
the same algorithm for the calculation of undisturbed sedi- 
ment temperatures (using the frictional decay) and thermal 
conductivity of the sediment (using the heat-pulse decay). 
Inversion theory allows calculation of realistic errors in a 
well-defined and mathematically rigorous way based on the 
sample rate and temperature resolution used. Contrary to 
HFRED, the thermal conductivity and diffusivity are treated 
as independent para meters, which may allo w improvement 
of the relationship of iHvndman et alj (1 19791) . 

In the following we propose an inversion scheme, 
coined T2C (Temperature to conductivity), that was used to 
process the heat flow measurements made on the German RV 
SONNE on the eastern flank o f Juan de Fuca during cruise 
SOI 1 1 in the summer of 1996 (Villinger & Fahrtteiln ehmeit 
119961) . In some concluding remarks we summarise our expe- 
rience up to now with the new reduction scheme. 



1 THEORETICAL BACKGROUND AND 
INVERSION SCHEME 

The theoretical background for the analysis of heat flow 
measurements described in the Introduction is discussed in 
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Bulla rl (11954. [Listed (1 19791) . iHvndman et all (Il979h . and 
Villinger & Davis d!987l)~ The following simplified model 



for the sensor string is used: A cylinder of radius a and infi- 
nite extent in the z-direction is situated in a homogenous in- 
finite material. Whereas the material surrounding the cylin- 
der has a finite thermal conductivity k and thermal diffusiv- 
ity k, the cylinder itself is of infinite conductivity and dif- 
fusivity, with the constraint that (pc) c , the product of spe- 
cific heat c and density p of the cylinder, remains finite. At 
time t — 0, the cylinder is at temperature Tq and the ambi- 
ent space at T a . The temperature at the centre r — of the 
cylinder can then be described by the thermal decay curve of 
the cylinder: 

T(T) = (T -T a )-F(a,r)+T a (1) 

with F(cx,t) follo wing iBullardl dl954 and 
ICarslaw & Jaeger! dl959l) : 



F(a,r) 



4a 

72 



du 



ix" J u<j)(u, a) 
o 

(uJo(u) — a,Ji(u)) 2 
+ (uY (u) - aYi(u)) 2 



(2) 



(3) 



Here a = 2(pc) s / (pc) c is the ratio of the heat capacities 
per volume of the sediment and the cylinder, and r = is 
the dimensionless time. Jo, J\, Yq and Y\ are zero and first 
order Bessel and Neumann functions, respectively and u is 
the integration variable. The temperature rise To — T a at time 
t = can be expressed by 



T Q - T a = 



Q 



na 2 (pc) c 



(4) 



with Q being the heat per unit length contained in the cylin- 
der. EquationQ]has an asymptotic solution th at approximates 
F(a t) with 1% accuracy for r > 10 dHvndman et all 
[1979) : 



T(t) 



Ankt 



(5) 



It is important to recall the limitations of this model by com- 
paring it with the sensor tube housing the temperature sen- 
sors: 

(i) The sensor tube is non-ideal, i.e. it has a finite conduc- 
tivity and an internal structure. 

(ii) The duration of the heating pulse is finite, usually in 
the order of 10-20s. 

(iii) Axial heat flow will be inevitable but certainly small. 

(iv) A thin water-layer between the sediment and the sen- 
sor tube may act as insulation to delay the achievement of 
thermal equilibrium. 

Measurements early in the temperature records will be in- 
herently affected by deviations from the model. Therefore, 
these records have to be excluded from analysis. However, 
temperatures within the analysed time range still show slight 
deviations from the ideal behaviour. This deviation can be 
best modelled by introducing a new parameter, the time shift 
t s . The measured time origin is always the onset of the pen- 
etration or heat-pulse. Introduction of the parameter t s ap- 
proximates the heating of finite length by an instantaneous 



temperature rise, shifted relative to the onset of the heat- 
ing. Although mathematically not rigorously proven, this 
concept is reasonable from a physical po int of view and 
has been shown to provide reli able results dHvndman et all 
1 19791: IVillinger & Davisl |1987|) . The justification for using 
t s will be investigated more thoroughly in section [4] using a 
numerical model. 

The goal of the processing scheme is to invert equa- 
tion Q] To achieve this, the actual physical parameters have 
to be restored in equations [2] and [3] This yields 



T 



-u 2 -^(t-t 3 ) 



ir 3 a 2 K(pc)c J u<j){u) 
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du + T n 



(6) 
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+ uY (u) 



2k 



n(pc). 



-Yi(u 



(7) 



The equation assumes that the initial temperature Tq at t = 
is reached by introducing an amount of heat Q either due 
to frictional heating during penetration or a calibrated heat 
pulse. 

According to equation (O the temperature decay is a 
function of six parameters and time: 

T = T{m,t) with m = [k,K,T a ,Q,(pc) c ,t B ] (8) 

The inversion will be demonstrated for the general case of 
all six parameters being unknown. In practical application of 
the technique, the number of parameters needs to be reduced. 
This will be discussed in section [2] 

The decay function is nonlinear and cannot be inverted 
directly. One approach to solve this problem is to expand 
T(m , t) in terms of a first order Taylor series (Menke, 
1989; Kristiansen, 1982). The temperature will be expanded 
around the "true" set of parameters r, using an estimated set 
of parameters e: 



T(r,t)=T(e,t) + J2 



dT 
drrii 




(9) 



This equation is linear in the difference vector A and can 
be inverted for this parameter vector. The left side of the 
equation represents the data and the right side the model. 
Because the equation is only of first order and hence not 
exact, a recursive scheme will be used for the calculation 
of the true parameters whereby the result of the Z-th iteration 
is used as an estimate e for the (I + l)-th iteration. 



e ('+D = e «> 



(10) 



In general e^ +1 ^ will be closer to the true parameter set 
than e(", Because the parameters are of different orders of 
magnitude, the difference vector A has to be normalised 
so that only relati ve weights appear in the model kernel 
dKristiansenl[l982l) . 



A' 



1. 



,6 



(11) 



Applying this method to all data points for times tj (j — 
1, . . . , N) yields a system of N linear equations. These can 
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Figure 3. Spectrum of the SVD for a representative set of parameters. For 
i > 4, magnitudes are in the range of floating point accuracy of the imple- 
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Table 2. Values of the parameters used to calculate a model kernel for fric- 
tional and heat pulse decay. Errors are calculated from equation [T9l For the 
fixed parameters in the inversion no error is given. 



be written in matrix notation: 



GA' 



dT(m,t s ) 



T(r,t 3 )-T(e,^ 



(12) 
(13) 
(14) 



G is the model matrix or kernel, it contains no actual data, 
only the assumptions of the model. A is the model parameter 
vector and d is the data vector. 

The solution of the problem requires the inversion of 
the matrix equation [T2] in order to obtain A. We use sin- 
gular v alue decomposition (SVD) closely following iMenkd 
(1989). The decomposition of the model matrix produces 
three matrices G = U • S • T T with the solution to the 
inverse problem: 



A n = VS" 



(15) 



Once conductivity- and temperature-depth profiles are 
calculated from the decay curves, heat flow values are com- 
puted using a method proposed bv lBullardl dl939l) . Starting 
from the heat diffusion equation for the horizontally layered 



case, 

q = ^-(k(z)T(z)) 
oz 

an integration over depth z yields: 

dz' 



T(z) - To 



k(z') 



(16) 



(17) 



The integral is called the thermal resistance w and is calcu- 
lated from the conductivity-depth profile. A linear regression 
of the temperature profile versus the thermal resistance gives 
the heat flow density as the slope of 



T(z) — q- w + T 



2 PROCESSING SEQUENCE 



(18) 



In the last section a general inversion scheme was proposed 
which is capable of inverting frictional and heat pulse de- 



cays. It is important that all parameters to be determined are 
independent of each other. A measure of the linear depen- 
dence is the magnitude of the diagonal values of the sin- 
gular matrix S. A plot of these values is called the spectrum 
of the singular value decomposition and is shown for a typ- 
ical temperature decay in figure [3] The accuracy of floating 

point figures given by Matlab® is about 10 -16 . The spec- 
trum drops to this level after the fourth singular value, sug- 
gesting that only four independent parameters can be deter- 
mined in the problem. However, the accuracy and stability of 
the inversion is largely influenced by the smallest included 
singular value and error bounds are proportion al to the in- 
verse of the square of this value (lMenke[ [l989). Therefore, 
although it is possible to use four independent parameters, in 
practice only three parameters were determined to improve 
the results. A thorough analysis revealed that the singular 
value A4 is mostly influenced by (pc) c . This observation is 
a result of the heat capacity being a probe parameter that is 
only poorly resolved by the model. Therefore a fixed value of 
3.38 -10 6 J/(Km 3 ) was used, base d on the material an d ge- 
ometry of the sensor string (Nagihara & Listerl fl993). The 
possibility of a systematic error resulting from a false as- 
sumption for the heat capacity will be discussed later. 

The restriction to three parameters for one inversion run 
leads to a c ertain sequence of proces sing steps for a single 
penetration dVillinger & DavisLl 19871) : 

(i) Frictional decay: Q,T a ,t s are calculated. k,K,(pc) c 
are held fixed, using estimated values for k and n. 

(ii) Heat pulse decay: The residual of the frictional decay 
is subtracted from the heat pulse decay thus reducing the am- 
bient temperature T a to a known value of zero. The inversion 
is used to compute k, k, t s with T a , (pc) c , Q held constant. 
Q is known as an input parameter and T a is zero after the 
reduction. 

(iii) The calculations in steps (i) and (ii) can be repeated, 
this time using the calculated values of k and n as improved 
estimates in step (i). 

The technique was first used on the results of research cruise 
SOI 11 (see section Penetrations took place over young 
crust and comparatively warm sediments. Especially for the 
upper thermistors on the sensor string it was noticed that the 
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Figure 4. Distribution of conductivity values when synthetic, normally dis- 
tributed noise with a = 1 mK is added to a synthetic data set. The mean 
value is (1.001 ± 0.017) W/(m K). 



frictional heat of the penetration did not suffice to raise the 
temperature from seawater temperature above the ambient 
sediment temperature. This caused very small signals for the 
frictional decay of some of the sensors, for example Til in 
figure |2] 

For the model we used it turns out that such a small fric- 
tional decay can be described equally well by a large heating 
shifted by a large amount in time or a small heating and a 
small time shift. This ambiguity in the parameters Q and t s 
cannot be resolved without additional information and leads 
to problems when inverting frictional decays. Because the 
important parameter T a is not influenced by this ambiguity, 
a practical approach to this problem is to set the time shift to 
a fixed value, namely zero for these cases. 

To determine occurences of this problem, the first or- 
der approximation (equation [5]l is used to determine decay 
curves with very little frictional heat. If Q, computed from 
this equation, is below a certain threshold, t s is set to zero. 

The stability of the inversion of heat pulse decays is af- 
fected by K-values that fluctuate from iteration to iteration. 
In the worst case this causes the algorithm to diverge. To pre- 
vent this, only small changes in n are allowed between iter- 
ations, effectively damping any oscillations. This constraint 
is used only for the first three iterations. 

3 CALCULATION OF ERRORS 

For inversions th eoretical erro rs can be calculated using the 
model kernel G jMenkelll989h : 

cov m = G _I cov d G- IT (19) 

Here cov m and covd are the covariance matrices of the 
model parameters and the data, respectively. G _I is the in- 
verse of the kernel, in our case VS _1 U T (equation [T5l>. 
This equation is only exact for linear problems with nor- 
mally distributed data. Nevertheless, it can be applied to a 
nonlinear problem, if the problem can be approximated by 
a linear function in the vicinity of the solution. It is a useful 
property of equation [19] that no data is used for the calcula- 
tion of errors, only the data covariance and a model kernel 
are needed. The kernel can be computed using equation [121 
by assuming a parameter vector with representative values. 



This choice is not critical since the kernel and correspond- 
ing errors vary only slowly with the parameters. This feature 
of equation [151 is particularly useful for design studies and 
theoretical evaluation of the limitation of a certain model. 

The covariance matrix of the data is not known a priori 
but in our case it can be assumed that data errors are uncor- 
related and mainly due to the finite temperature resolution of 
±1 mK. covd is then a matrix with temperature variances on 
its main diagonal and zeros elsewhere. To compute the ker- 
nel G using equation[12] a modelled decay curve from 120- 
420 s with 10 s sampling rate and two sets of parameters for 
frictional and heat pulse decay are used, respectively. These 
values, together with the calculated errors, are summarised 
in table [2] For the frictional decay, errors for Q,t s and T a 
are determined, according to the processing sequence given 
in section [2] The ambient temperature can be resolved with 
a standard deviation of 1.4 mK, or with a relative error of 
0.28%. For the heat pulse decay, errors for k, K and t s were 
computed. The standard deviation of the thermal conductiv- 
ity is 0.018 W/(mK) (1.8%). 

The error bounds given in table [2] give a general idea 
of the errors to be expected when using the algorithm. The 
relevant parameters conductivity and sediment temperature 
can be calculated with relative errors of about 2.0 and 0.5%, 
respectively. For the parameter range encountered in prac- 
tice, errors will vary slightly. Therefore, for each penetration 
errors are calculated from the inversion results and stored 
together with them. 



4 TESTS OF THE PROPOSED INVERSION 
SCHEME 

To verify the theoretically derived error bounds, normally 
distributed noise with a = 1 mK was added to the exact 
model data. This data was then inverted. To obtain rea- 
sonable statistics the procedure was repeated 100 times, 
producing results for k which are shown in the histogram 
in Figure [4] The mean of the computed conductivities, 
(1.001±0.017) W/(mK), agrees very well with the ex- 
pected value of (1.00±0.018) W/(m K). 

The value of (pc) c is only known on the b asis of geo- 
metric considerations (Nagihara & Listed . Il993l) . Addition- 
ally, its value could vary by a few percent over the length of 
the sensor string because of more wiring in the upper parts 
of the string. If a wrong value for (pc) c is assumed it will 
influence the results of the inversion in a systematic way. 
Therefore it is useful to investigate the magnitude of the er- 
ror introduced by an incorrect value. 

A synthetic temperature record was calculated using 
the same parameters as in the previous section. (pc) c 
was chosen to be 3.38 T O 6 J/(Km 3 ), the value given by 
Nagi hara & Listed dl993l) . This record was the input to an 
inversion in the heat pulse-configuration (fc, k, t s are cal- 
culated). Several runs were made, each time with a slightly 
different value for the heat-capacity, ranging from 2.9 to 
3.9T0 6 J/(Km 3 )which corresponds to a relative deviation 
of about ±15%. 

Figure [5] shows relative errors in the calculated model 
parameters when (pc) c is varied. For the maximum bias, the 
deviation of the conductivity is only about 0.2%. This value 
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Figure 5. Influence of false (pc) c -values on inversion results. Variations of 
up to 15% result in no significant change in k. 
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Table 3. Inversion results of the TFELD model. The results show good 
agreement. 




Diffusion time 



is considerably smaller than the errors introduced from other 
sources. 

The situation is different for the other two parameters 
though. The systematic deviations are significant compared 
to the random errors. t s only has meaning within the frame- 
work of the inversion algorithm, but n is a physical parame- 
ter and a relationship based on the calculated k could be in 
error by a few percent. 

The tests were conducted using synthetic data based on 
equation [6] A time shift can be computed for model curves 
but it cannot be verified whether t s is able to approximate 
the non-ideal parts of real temperature records. 

A closer approximation to reality is a numeric model 
that is able to model finite heat pulse length and finite probe 
conductivity. For this reason a program called TF ELD was 
used w hich is based on an algorithm proposed bv lVillirigerl 
(fl985h . This program allows modelling of heat diffusion in 
a cylindrically layered space with various heating functions. 
Each of the cylindrical layers is characterised by a set of 
physical parameters p, k and c. In this case a model with 
only two layers was employed, the inner and the outer layer 
representing the probe and the sediment, respectively. 

In the first model used the conductivity and diffusiv- 
ity of the probe (layer 1) were set to 1000 W/(mK) and 
2000- 10 -7 m 2 /s respectively to give a good approximation 
of an ideal probe. The ambient sediment (layer 2) values 
of 1 W/(mK) and 2-10~ 7 m 2 /s were used as representa- 
tive values for conductivity and diffusivity. An instantaneous 
temperature rise of the inner cylinder was used as the heat 
source. The computed temperature decays were used as in- 
put for the inversion algorithm. The results are shown in 
table [3] As expected the inversion gives the correct result 
within the precision of the numerical computations. 

The main objective of using the TFELD model was to 
study the influence of a non-ideal probe on the tempera- 
ture decay curves . In the literature dHvndman et all 119791: 
IVillinger & DavisL 119871; iNagihara & Listen, 1 1993b it is as- 
sumed that the effects of the finite thermal conductivity of 
the cylinder and the finite heat pulse on the temperature de- 
cay will be combined in the time shift t s . 

In order to separate effects of finite heating and probe 
geometry two experiments were conducted, one with an ap- 
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Figure 6. Influence of non ideal probe parameters: a) finite heating b) finite 
diffusivity. 



proximately ideal probe and finite heat pulse and a second 
with real probe geometry and infinitely short heat pulse. 

For the first experiment probe conductivity was fixed at 
1000 W/(m K) and heat pulses with duration from 10-20 s 
were used as input to TFELD. The resulting decay curves 
were inverted using T2C. If the time shift from the inversion 
is plotted versus the length of the heat pulse (Figure |6^), a 
strong linear relationship can be seen. 

A similar test was run to determine the relationship be- 
tween probe geometry and time shift. It was found that the 
parameter suited best to describe the probe geometry is the 
diffusion constant td'. 

a 2 

t d = — (20) 

K 

This value describes the time that a temperature disturbance 
would take to travel the distance a. The diffusion time of 
the probe was varied from 2^4-0 s in this experiment cor- 
responding to values for k and k of 40-2 W/(mK) and 
80-4 -10 -7 m 2 /s, respectively. The initial temperature field 
was To inside the probe and zero outside. Again these model 
curves were inverted and the resulting time shifts plotted ver- 
sus the diffusion time to obtain the linear relationship seen 
in Figure |6p. 

If a linear fit is calculated for both curves, the following 
relationships are obtained: 

t a = 0.554/, - 0.56 [s] (21) 
t a = 0.3lt d - 0.08 [s] (22) 

Here th is the duration of the heat pulse. It is intuitively ex- 
plicable that a temperature record, generated by a finite heat 
pulse can be best described by a perfect pulse shifted half 
the length of the original pulse. Similarly the probe geome- 
try will result in a time shift that is directly dependent on the 
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time the heat pulse needs to travel through the internal struc- 
ture. This is reflected in equations [21] and |22j respectively. 

One reason for the development of our algorithm was to 
evaluate the possibilities of calculating the conductivity di- 
rectly from the frictional decay. In this case the only known 
parameter is (pc) c . From the discussion of the spectrum of 
the S VD, it is clear that a reduction of the number of param- 
eters to at least four is necessary to invert equation Q] as an 
overdetermined problem. To reduce the number of parame- 
ters in the problem by one, a relat ionship between co nduc- 
tivity and diffusivity can be used dHvndman et all [l 9791) . 



Thermal parameters: 



5.79 - 3.67fc+1.016/c 2 



10" 



(23) 



This assumption is justified because both conductivity and 
diffusivity of marine sediments are mainly porosity con- 
trolled. Using equation [19] one can calculate expected error 
bounds for this configuration. If an observational error of 
10~ 3 K and reasonable parameters for probe and sediment 
are assumed, the errors in the model parameters are as fol- 
lows: 

A/c = 6.6W/(m K) 
AQ = 1500 J/ m 
At s = 58 s 
AT a = 12 mK 

It can be seen in equation [19] that the model covariance is 
linearly related to the data covariance. If an error margin of 
0.05 W/(m K) (w ±5 % for marine sediments) is assumed 
as an acceptable error level for the conductivity, this means 
a reduction of Afc by 2 orders of magnitude. This means 
as well a reduction of the observational error by 2 orders 
of magnitude. Thus, accuracy of temperature readings has 
to be reduced to an extremely low and unrealistic level of 
±10~ 5 K. The results of these tests confirm that the fric- 
tional decay is not sensitive to the thermal conductivity of the 
surrounding sediment and a successful inversion for k with 
realistic error bounds for the temperature measurements is 
not possible. 

The reason for this failure is implicit in the model. Not 
enough information is provided by a frictional decay curve 
to resolve the ambiguity of the problem using an inver- 
sion scheme. Programs exist, however, that are capable of 
modelling the full fricti onal decay (forward modelling), in- 
cludin g the early times dLee & Von Herzenlll994tlwi inger. 
1 1985b . This allows other constraints to b e added to the prob- 
lem, thereby reducing the ambiguity dLee & Von HerzenL 
1 1 9941) . It has to be considered, however, that any additional 
assumptions and a-priori information will directly affect 
computation of thermal conductivities and have to be cho- 
sen carefully. 



5 INVERSION OF NEEDLE-PROBE 
MEASUREMENTS 

To test the algorithm with real data, a set of pulsed needle- 
probe measurements were used. These were made in a ce- 
ramic alloy with a known thermal conductivity close to the 
value of deep sea sediments. The data were kindly provided 
by TeKa Inc. (Berlin, Germany), a company specialised in 
thermal conductivity measurements and equipment. The true 



Conductivity (1.609 ± 0.016) W/(m K) 



Diffusivity 



7.3-10" 7 m 2 /s 



Heat capacity 2.20- 10 6 J/(K m 3 ) 



Measurement parameters: 



Heating power 8, 12, 15 W/ m 



Heat pulse duration 5, 10 s 



Sampling rate 0.5 s 



Observational error 10" 4 K 

Table 4. Upper half: Thermal parameters of the material used for needle- 
probe measurements. Lower half: Parameters of the different records taken. 



# ^ # ^ ^ ^ 



6 6 



i 



1.625 
1.62 
| 1.615 
— 1.61 
1.605 
1.6 



Figure 7. Inversion results for the needle-probe measurements. The circles 
represent calculated conductivities with 1-<t error bars for the 18 records. 



conductivity of the ceramics was measured with a divided 
bar device (TeKa Inc., pers. comm.). The relevant param- 
eters of the material and the recording settings are sum- 
marised in table[4] For each combination of parameters three 
decay curves were recorded, resulting in a total of 18 decay 
curves. 

A drift correction was applied to the temperature 
records to reduce the value of the ambient temperature T a 
to zero. The heat capacity of the probe was unknown so that 
the decay curves were inverted using four unknown param- 
eters: [k, k, (pc) c , t s }. The results for the conductivity with 
1-ct error bars are shown in figure [7J The circles represent 
the computed thermal conductivities together with 1-<t er- 
ror bars. The dashed lines are the mean thermal conductiv- 
ity and 1% error margins determined with the divided bar 
method. Shown are all 18 measurements in six groups of 
three measurements with the same heating parameters. On 
the top horizontal axis, the heating parameters are encoded 
in the names: The first two digits represent the duration [s], 
the last two digits the heating power per unit length [W/m] 
of the heat pulse. For each set of parameters, three measure- 
ments were made. 

There is excellent agreement between expected and cal- 
culated values of the heat conductivity. Errors of k were be- 
tween 0.2 and 0.5%. For the K-values the errors are 3-15% 
and for {pc) c 10-30%. The high accuracy of the conductivity 
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0.6 
0.5 



Thermal Conductivity [W/(m K)] 
0.8 1 1.2 1.4 1. 



I 2 



8.5 9 9.5 10 10.5 11 11.5 12 12.5 13 
t s [s] 

Figure 8. Distribution of calculated time shifts for the needle-probe mea- 
surements. Two clearly distinguishable groups can be identified that are 
about 2.5 s apart from each other. 



value compared to previous tests is due to the small obser- 
vational error of the needle probe. As stated before, large 
errors for k and (pc) c are related to the inversion with four 
unknown parameters. 

The error bounds generally decrease to the right in fig- 
ure [7] Changes in length of the heat pulse and heating power 
affect the total heat input and the recorded temperature sig- 
nal. Thus it appears that calculated error bounds vary with 
the magnitude of the recorded temperature changes, as ex- 
pected. 

It is interesting to check whether the hypothesis that the 
time shift is connected to the diffusion time via a simple re- 
lationship of the type in equation[2T|holds true for real data. 
Figure [8] shows that all t s fall into two clearly distinguish- 
able groups according to the duration of the heat pulse (5 and 
10 s). The difference of about 2.5-3.0 s is in good agreement 
with the expected value, which should be half the difference 
of the lengths of the heat pulses. 



6 RESULTS FROM RESEARCH CRUISE SOlll 

During summer of 1996 cruise SOlll was conducted off 
Western Canada on the German research vessel SONNE. 
The objective was to study the effect of hydrothermal cir- 
culation on marine heat flow on the Eastern flank of the Juan 
de Fuca Ridge. The survey area was located near the Cobb 
Offset at about 47° 30' N, 129° 0' W. During the cruise 8 sta- 
tions with 104 successful penetrations were made. For a de- 
tailed discussion of the measureme nts during this cruise see 
IVillinger & Fahrtteilnehmerl d 19961) . 

This cruise provided the first instance to use the pro- 
gram on a regular basis. The following is concerned mainly 
with the application of the inversion algorithm. Though all 
penetrations were inverted successfully using the described 
algorithm only data from station 2 with 30 penetrations will 
be used to illustrate the applicability of the program. 

In the first stage of processing the moment of penetra- 
tion of the probe and the onset of the heat pulse have to be 
picked manually using the raw temperature data. The data 
can then be input to the inversion program that calculates 
thermal parameters according to the previously described se- 




Figure 9. Conductivity-depth profile for station 2. 



quence. The thermal parameters are then in turn used to cal- 
culate heat flow values. 

To minimise the influence of the non-ideal probe at early 
times, only temperatures for times > 120 s relative to the 
picked origin time were used. The value for this start time 
depends on the probe geometry and was found empirically 
by examining the fit for several start times. For example 
Nagi hara & Lister! (Il993l) used a sensor string with a diame- 
ter of 9.52 mm versus 8.0 mm in our design and a value of 
200 s as the starting time for their analysis. 

After calculation of conductivities and thermal gradient, 
the absolute penetration depth can be calculated using the 
bottom water temperature and the assumption that the tem- 
perature is continuous at the sediment-water interface. Fig- 
ure|9]shows the conductivity versus the sensor depth for sta- 
tion 2. A strong increase in conductivity over the first metre 
can be seen that is possibly caused by a decrease of poros- 
ity in the surface sediments. The outliers in the lower part 
of the plot might be attributed to a partially penetrated sand 
layer of varying depth. This corresponds to the observation 
that often the penetration was stopped by a layer of harder 
material. 

Previous workers (iBullardl (|1954|): iRatcliffd (T l 960 ) ; 
IVon Herzen & Maxwell ldl959h : lHvndman et al.ldl979l) ) de- 
termined thermal diffusivity by measuring conductivity, den- 
sity and porosity. Assuming that heat capacity and conduc- 
tivity are controlled mainly by the water content, they con- 
structed a relationship between thermal conductivity and dif- 
fusivity. In our approach, thermal diffusivity is computed as 
an independent parameter in the inversion. Figure [TOl shows 
calculated diffusivity v alues compared t o the relationship 
(Equation |23l given by iHvndman et al] d!979l) . The inver- 
sion results are significantly lower than expected by the the- 
oretical curve. To our knowledge no other measurements of 
our type for marine sediments are published to verify the 
results. It is apparent, though, that a simple equation is not 
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k [W/(m K)] 

Figure 10. Calculated diffusivity versus conducti vity for station 2. The v al- 
ues are compared to an empirical relationship by Hyndman et al. 1 1979). 
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k(HFRED) [W/(m K)] 



Figure 11. Comparison of thermal conduct ivities between T2C described 
in this paper and HFRED IVillinger & DavisUl987l) . Dashed lines show a 
relative deviation of ± 5%. 



sufficient to describe the variations in the relation between k 
and k. 

It is instructive to compare the results of HFRED 
dVillinger & Davisl 119871) to our results. Figure [TT| shows 
differences in computed conductivity for both algorithms. 
A high agreement is visible with deviations usually within 

±5%. 

HFRED uses a constant value for the ratio of the heat 
capacities a of 2.0 because of the lack of knowledge of the 
thermal diffusivity. In the T2C algorithm a is not used but 
can be calculated from the computed thermal diffusivity as 
k = k/(pc) c . Figure [12] shows relative differences of the 
thermal conductivity versus a. It is apparent that large devi- 
ations from a = 2 correspond to large differences in thermal 
conductivity, suggesting that a constant value of a is not suf- 
ficient to describe all measurements accurately. 



10 




a 



Figure 12. Relative differences of conductivities computed with HFRED 
and T2C as a function of a. 



7 CONCLUSION 

An algorithm is described to invert heat flow measurements 
made with a violin-type heat probe including in-situ ther- 
mal conducti vity measurements after the pulse method of 
iListerl (119701) . A theoretical analysis of the inversion al- 
gorithm shows that undisturbed sediment temperatures can 
be determined with an error of about 1-2 mK. The abso- 
lute error of the in situ thermal conductivity k is less than 
0.002 W/(m K). It is also possible to compute thermal dif- 
fusivities, but with considerably higher errors of about 10 %. 

Numerical experiments with synthetic temperature de- 
cay data reveal a strong relationship between the time shift 
and the thermal probe parameters which can be explained in 
a quantitative way by the finite length of the heat pulse and 
the diffusion time constant of the sensor string. 

Measured data, obtained with a needle probe in mate- 
rial with known thermal conductivity, confirm the accuracy 
of the inversion procedure and show that the algorithm is 
suited to the analysis of pulsed needle probe measurements. 
Our results show that is is possible to succesfully use the 
algorithm in pulsed needle probe measurements. 

The inversion of the data obtained on SOI 1 1 proves that 
the described algorithm is robust and well suited for auto- 
mated processing of a large number of heat flow penetra- 
tions. The embedding of the software in a suite of mathe- 
matical software allows simple further analysis of the data 
and easy development of additional tools. The relative ac- 
curacy of our thermal conductivity results is in the range 
of 1-3%. Undisturbed sediment temperatures can be com- 
puted with relative errors of 0.5-1%. The comparison with 
results obtained wi t h the previously used program HFRED 
dVillinger & Davisl 1 1987b shows good agreement between 
both algorithms. Deviations are generally due to the assump- 
tion of a = 2 by HFRED. 



References 

Bullard, E. C, 1939. Heat flow in South Africa, Proc. R. 

Soc. London, Series A, 173, 474-502. 
Bullard, E. C, 1954. The flow of heat through the floor of 

the atlantic ocean, Proc. R. Soc. London, Series A, 222, 

408^125. 



10 A. Hartmann and H. Villinger 



Carslaw, H. S. & Jaeger, J. C, 1959. Conduction of heat in 
solids, Oxford University Press, 2nd edn. 

Hyndman, R. D., Davis, E. E., & Wright, J. A., 1979. The 
measurement of marine geothermal heat flow by a multi- 
penetration probe with digital acoustic telemetry and in 
situ thermal conductivity, Mar. Geophys. Res., 4, 181— 
205. 

Jones, E. J. W., 1999. Marine Geophysics, Wiley, Chich- 
ester, UK. 

Kristiansen, J. I., 1982. The transient cylindrical probe 
method for determination of thermal parameters of earth 
materials, GeoSkrifter 18, Laboratory of Geophysics, 
Aarhus, Denmark. 

Lee, T.-C. & Von Herzen, R. R, 1994. In situ determination 
of thermal properties of sediments using a friction-heated 
probe source, J. Geophys. Res., 99(B6), 12121-12132. 

Lister, C. R. B., 1970. Measurement of in situ sediment 
conductivity by means of a Bullard-type probe, Geophys. 
J. R. astr. Soc, 19, 521-532. 

Lister, C. R. B., 1979. The pulse-probe method of conduc- 
tivity measurement, Geophys. J. R. astr. Soc, 57, 451- 
461. 

Menke, W., 1989. Geophysical Data Analysis: Discrete In- 
verse Theory, no. 45 in International Geophysics Series, 
Academic Press, rev. edn. 

Nagihara, S. & Lister, C. R. B., 1993. Accuracy of marine 
heat-flow instrumentation: Numerical studies on the ef- 
fects of probe construction and the data reduction scheme, 
Geophys. J. Int., 112, 161-177. 

Ratcliffe, E. H., 1960. The Thermal Conductivity of Ocean 
Sediments, J. Geophys. Res., 65(5). 

Villinger, H., 1985. Solving cylindrical geothermal prob- 
lems using the Gaver-Stehfest inverse Laplace transform, 
Geophysics, 50(10), 1581-1587. 

Villinger, H. & Davis, E. E., 1987. A New Reduction Al- 
gorithm for Marine Heat Flow Measurements, J. Geophys. 
Res., 92(B12), 12846-12856. 

Villinger, H. & Fahrtteilnehmer, 1996. Fahrtbericht SOI 1 1, 
Berichte aus den Geowissenschaften 97, Universitat Bre- 
men. 

Von Herzen, R. & Maxwell, A. E., 1959. The Measure- 
ment of Thermal Conductivity of Deep-Sea Sediments by 
a Needle-Probe Method, J. Geophys. Res., 64(10), 1557- 
1563. 



ACKNOWLEDGMENTS 

The paper was improved by the thorough reviews of our 
colleagues, Ingo Grevemeyer, Norbert Kaul, and Marion 
Pfender as well as R. von Herzen and an anonymous re- 
viewer. The research cruise SOI 11 was kindly funded by 
the German Ministry for Research and Technology, Grant 
No. 03G0111A. 



